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Abstract 

A mechanism of Mott transitions in a Bose Hubbard model on a square lattice is studied, using a variational 
Monte Carlo method. Besides an onsite correlation factor, we introduce a four-body doublon-holon factor into 
the trial state, which considerably improves the variational energy and can appropriately describe a superfluid- 
insulator transition. Its essense consists in binding (and unbinding) of a doublon to a holon in a finite short 
range, identical with the cases of fermions. The features of this transition are qualitatively different from those of 
Brinkman-Rice-type transitions. 
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1. Introduction: After early theoretical studies 
of Mott or superfluid-insulator transitions in inter- 
acting Bose systems [1] , an experimental example 
has been realized using an ultracold dilute gas of 
bosonic atoms in an optical lattice [2] . The essence 
of this system is considered to be captured [3] by a 
Bose Hubbard model. This basic model has been 
studied with various methods; for square lattices, 
properties of T c , superfluid density, etc. were stud- 
ied, applying a quantum Monte Carlo method to 
small systems (mainly 6x6 square lattice) [4], 
and a ground-state phase diagram in a plane of 
chemical potential and interaction strength was 
obtained, using a strong-coupling expansion [5]. 
These studies estimated the critical interaction 
strength of Mott transitions at U c /t = 16.4 ± 0.8 
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and 16.69, respectively, for the particle density 
of n = 1 (n = N c /N with N c : particle number, 
N: site number) at T = 0. Thus, the existence 
of a Mott transition has been embodied, but the 
mechanism of the transition is still not clear. 

Variational Monte Carlo approaches are very 
useful to understand the mechanism of the Mott 
transition, because one can directly and exactly 
treat wave functions. For the Bose Hubbard model, 
wave functions with only onsite correlation factors, 
which corresponds to the celebrated Gutzwillcr 
wave function (GWF, Wg) for fermions [6], were 
studied first [7]. In contrast to for fermions, GWF 
for bosons is solved analytically without additional 
mean-ficld-type approximations [8] for arbitrary 
dimensions, and yield a Brinkman- Rice- type (BR) 
transition [9] at U = Ubr- In BR transitions, all 
the lattice sites are occupied with exactly one 
particle and the hopping completely ceases in the 
insulating regime, namely, ^g ~~ * IljLi bj\0) an d 
E = for U > Ubr- This result is caused by 
an oversimplified setup of the wave function, in 
which the effect of density fluctuation should be 
included. In this work, we introduce a doublon- 
holon binding correlation factor into the trial 
function, following previous studies for fermions. 
Thereby, we can describe a superfluid-insulator 
transition more appropriately. 

2. Formulation: We consider a spinless Bose- 
Hubbard model on a square lattice, 

H = -tJ2 (b\bj + b]h) + ^ ]T - !)> C 1 ) 

(*?'} 3 

where bj is a creation operator of a boson at site 

j, rij = b^bj and t, U > 0. Here, (ij) denotes a 
nearest-neighbor-site pair; the definition of t is a 
half of what was given in some literatures [4,7]. In 
this work, we restrict ton~ 1, because we would 
like to consider the most simple case of Mott tran- 
sitions. The cases of other commensurate densities 
(n > 2) must be essentially identical. 

We study this model eq. (1) through a varia- 
tion theory. As a trial wave function, we use a 
Jastrow type of *q = PqP g $ (QWF), follow- 
ing previous studies for fermions [10,11]. Here, <3? 
is the ground state of noninteracting (completely 
coherent) bosons, namely <& = 1, Pq is an onsite 



projector corresponding to the famous Gutzwillcr 
factor for fermions [6]: Pg(<7) — 9 D > with D = 
J2i n i{ n i ~ l)/2- As we repeatedly showed [11,12], 
intersite correlation factors are indispensable for 
appropriate descriptions of interacting systems. In 
particular near half filling, a four-body doublon- 
holon correlation factor Pq is crucial to explain 
the mechanism of Mott transitions [13,14,15]. For 
S = 1/2 fermions, Pq is explicitly written as, 

P Q (p) = (l-^ (2) 
with < \i < 1 and 

Q = Y,I[[ d i( 1 - e i+r) + ei(l-d i+T )], (3) 

i t 

where di (= n^nn) and [= (1 — nj|)(l — 
arc the doublon and holon operators respectively, 
i runs over all the sites, and r the four nearest- 
neighbor sites of the site i. When the variational 
parameter fi vanishes, is reduced to the GWF, 
$G = Pg&, in which doublons and holons can 
move independently. On the other hand, in the 
limit of fi — > 1, a doublon becomes bound to a 
holon in nearest-neighbor sites. Consequently, plus 
(holons) and minus (doublons) charge carriers 
completely cancels, indicating a metal-insulator 
transition. 

In this work, we extend eq. (2) to bosons, and 
include the correlation between diagonal (second- 
nearest)-neighbor sites: 

Pq(M,m') = (1-m) Q (1-m') Q ', (4) 

where the primes (') denote the cases of diagonal 
neighbors. In eq. (4), Q has the same form as Q 
in eq. (3), but the doublon operator di is replaced 
by a multiplon operator nij, which yields 1 (0) if 
the site i is multiply occupied (otherwise). Near 
Mott transitions, Pq of eq. (4) is substantially a 
doublon-holon factor, because the probability den- 
sity of multiplon with more than two bosons is al- 
most zero for such large values oiU/t. 

To estimate expectation values accurately, we 
use a variational Monte Carlo method [16,17,18] of 
fixed particle numbers. First, we optimize the var- 
iational parameters, g, fj, and //, simultaneously 
for each set of U/t, n and L, and then calculate 
physical quantities with the optimized parameters. 
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Fig. 1. Comparison of total energy per site between GWF 
and QWF as a function of U It for several system sizes. 
The particle density is commensurate (n = 1.0). For GWF, 
a Brinkman-Rice-type transition occurs at U = Ubh- The 
dash-dotted line is a curve proportional to —t/U. 

Through the optimization process, we average sub- 
stantially several million samples, which reduces 
statistical errors in the total energy typically to 
the order of 10~ 4 i. To check system-size depen- 
dence, particularly near phase transitions, we em- 
ploy square lattices of N = L x L sites up to TV = 
1024 (L — 32) with the periodic-periodic bound- 
ary conditions. 

3. Results: We start with comparison of the total 
energies E/t between GWF and QWF, which are 
shown in Fig. 1 for n = 1. As mentioned, in GWF 
[7], a BR transition occurs at Usn/t = 12 + 8\/2 = 
23.31...; for U > Ubr, each site is occupied by ex- 
actly one particle and hopping or density fluctu- 
ation completely ceases. Consequently, E/t van- 
ishes in the insulating regime. On the side of Bose 
fluid (U < Ubr), E/t vanishes as oc (1 — U /Ubr) 2 , 
meaning this transition is a continuous type. The 
difference of the two functions is small for small 
U/t ( ^ 15), but it becomes conspicuous as U/t 
approaches Mott critical values (U/t ^ 20). Thus, 
QWF is a considerably improved function in the 
point of the variation principle. 

As we will discuss in detail shortly, QWF also ex- 
hibits a superfluid-insulator transition, but its be- 
havior is qualitatively different from that of GWF. 
In contrast to E of GWF, E of QWF in the insulat- 
ing regime (U > U c ) does not vanish but is propor- 
tional to —t 2 /U as seen in Fig. 1, which behavior is 
expected from strong-correlation theories, namely, 
the density fluctuation does not cease but is re- 
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Fig. 2. Magnification of total energy of QWF near Mott 
transitions arising at U = U c indicated by arrows. Data 
for five system sizes are plotted; for L > 20, double-mini- 
mum or cusp behavior is observed at U = Uc- Solid sym- 
bols denote the energies of the optimized states, and open 
symbols near U c metastable states. 
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Fig. 3. Behavior of doublon-holon binding parameter near 
U c - Data for five system sizes arc plotted; for L > 20, 
discontinuities are observed at U = U c - Solid and open 
symbols have the same meaning as in Fig. 2. 

stricted to a finite short range. This is the essence 
of the mechanism of Mott transition owing to the 
doublon-holon binding. In Fig. 2, we show the mag- 
nification of E/t of QWF for U <~ U c . For large 
L (> 20), we find, near U = U c , double-minimum 
structure of E/t in the space of variational parame- 
ters, meaning this transition is first order, although 
such structure cannot be confirmed for L < 14 by 
our VMC calculations. 

The first-order features are more easily found 
in the variational parameters and in some quanti- 
ties. In Fig. 3, we plot the optimized doublon-holon 
binding parameter /i near U c . For L > 20, there 
are clear discontinuities at U — U c ; the large val- 
ues of \x for U > U c indicate that a doublon and a 
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Fig. 4. Behavior of (D), which is substantially the doublon 
density and an order parameter of Mott transitions. Data 
are shown only for the optimized states. 
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Fig. 5. The occupation rate of the k = (0, 0) level for the 
two wave functions is depicted as a function of U /t. Data 
for several system sizes are simultaneously plotted. 

holon arc tightly bound in nearest-neighbor sites. 
In Fig. 4, we show the average of D, which is, near 
U c , virtually identical with the doublon density, 
namely, an order parameter of Mott transitions. 
The discontinuities of this quantity at U c corrobo- 
rate a first-order Mott transition. 

Finally, let us consider a couple of properties of 
this Mott transition. In Fig. 5, we show the occu- 
pation rate, p(k) = (b k bk) /N, of the lowest-energy 
level, k = = (0, 0) versus U/t. Although p(0) 
does not directly indicate the superfluid density, it 
must be a good index of superfluidity. For the non- 
interacting case (U/t = 0), all the particles fall in 
the k = level. As the interaction becomes strong, 
p(0) decreases at first gradually and drops discon- 
tinuously to the order of 1/N at U c . Thus, the su- 
perfluidity vanishes at the transition. In Fig. 5, we 
plot the chemical potential, £ = dE/dn, estimated 
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Fig. 6. The Behavior of chemical potential, £/t, as a func- 
tion of particle density near n = 1 is shown for two values 
of U/t, 16 in the superfluid regime and 24 in the insulat- 
ing regime. Data of three system sizes are simultaneously 
plotted for each U /t. The singular behavior in ( at n = 1 
for U > U c stems from the singularly low E at n = 1 as a 
function of n. 

from finite differences, as a function of n. In the su- 
perfluid regime slightly below U c (U/t = 16, gray 
symbols) , £ is a smooth function of n even at n = 
1, indicating the state is gapless in density excita- 
tion. On the other hand in the insulating regime 
slightly above U c (U/t = 24, black symbols), ( has 
a large discontinuity at n — 1 ; a density excitation 
gap opens for U > U c at n = 1. The gap behavior 
is also confirmed by the density correlation func- 
tion N(q) for small |q| (not shown). 

4- Discussions: In this proceedings, we have 
found that a wave function with doublon-holon 
correlation factor, ^q, qualitatively improve the 
description of a Mott transition also in a Bosc 
system. Thus, it is probable that the mechanism 
of Mott transitions for bosons is basically identi- 
cal to that for fermions. It is urgent to compare 
theoretical results with experiments particularly 
of optical lattices. We have left many issues to be 
discussed, which will be published elsewhere soon. 
When main calculations here were finished, we 
became aware that a similar wave function had 
been studied recently [19]. 
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